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Transport phenomena still stand as one of the most challenging problems in computational physics. 

By exploiting the analogies between Dirac and lattice Boltzmann equations, we develop a quantum 
simulator based on pseudospin-boson quantum systems, which is suitable for encoding fluid dynamics 
transport phenomena within a lattice kinetic formalism. It is shown that both the streaming and 
collision processes of lattice Boltzmann dynamics can be implemented with controlled quantum 
operations, using a heralded quantum protocol to encode non-unitary scattering processes. The 
proposed simulator is amenable to realization in controlled quantum platforms, such as ion-trap 
quantum computers or circuit quantum electrodynamics processors. 
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Transport phenomena in fluid flows play a crucial role for many applications in science and engineering. Indeed, a 
large variety of natural and industrial processes depend critically on the transport of mass, momentum and energy 
of chemical species by means of fluid flows across material media of assorted nature PJ. The numerical simulation 
of such transport phenomena still presents a major challenge to modern computational fluid dynamics. Among the 
reasons for this complexity stand out the presence of strong heterogeneities and huge scale separation in the basic 
mechanisms, namely advection, diffusion and chemical reactions Eli- In the last two decades, a novel concept for the 
solution of transport phenomena in fluid flows has emerged in the form of a minimal lattice Boltzmann (LB) kinetic 
equation. This approach is based on the statistical viewpoint typical of kinetic theory mu- LB is currently used 
across a broad range of problems in fluid dynamics, from fully developed turbulence in complex geometries to micro 
and nanofluidics IS 0, all the way down to lattice gas automata 0 and quark-gluon applications 0. 

Recent improvements in ion trap and superconducting circuit experiments make these platforms ideal for challenging 
quantum information and simulation tasks. Trapped-ion experiments have demonstrated quantum information and 
simulation capabilities )10H12j , including the quantum simulation of highly correlated fermionic systems m , fermionic- 
bosonic models 531 QS| and lattice gauge theories [16j . Superconducting circuit setups can host nowadays top-end 
quantum information protocols, such as quantum teleportation m and topological phase transitions [TS]. These 
quantum devices are approaching the complexity required to simulate both classical and quantum nontrivial problems, 
as proposed by Feynman some decades ago |19j . Efforts in designing quantum algorithms for the implementation of 
fluid dynamics make use of quantum computer networks EDI El]. In these works, the quantum degrees of freedom 
are used on the same ground as classical parameters, and the exponential gain of quantum computers is not properly 
exploited. In contrast, systems described by pseudospins coupled to bosonic modes, such as the aforementioned ion- 
trap and superconducting circuit platforms, can enjoy quantum superposition and have advantages with respect to 
pure-qubit quantum computers in simulating fluids. 

In this article, we propose a quantum simulation of lattice Boltzmann dynamics, using coupled pseudospin-boson 
quantum platforms. Based on previously established analogies between Dirac and LB equations, we define here a full 
quantum mapping of transport equations in fluid flows. The LB dynamics is simulated sequentially by performing 
particle streaming and collision steps. The non-unitary collision process can be implemented with an heralded protocol, 
by sequential collapses of an ancillary qubit. The proposed mapping is amenable to realization in trapped-ion and 
superconducting circuit platforms. 


RESULTS 

The lattice Boltzmann equation is a minimally discretized version of the original Boltzmann’s kinetic equation, in 
which the fluid is modeled as an ensemble of particles that move and collide within a uniform lattice. The lattice 
Boltzmann dynamics is described by the equation 

(d t +v*V b )fi(x,t) = -Aijlfj&t) - f- q {x,t)]. (1) 

Here, fi(x, t) is the ith component of the particle fluid density associated with the lattice site x at the time t, and with 
discrete velocity v). The macroscopic fluid density at the site x is retrieved as p(x,t) = JA while the fluid 

velocity is defined as the weighted sum of the discrete velocities, u(x,t) = 1/p JA fi(x, t)vi. The velocity components 
fiVi, with i = 1,2,. ..Q, satisfy mass-momentum-energy conservation laws and rotational symmetry. Typical lattices 
are D2Q9 or D3Q15 models, for the case of two dimensions with 9 speeds, and three dimensions with 15 speeds, 
respectively fTI \ . 

Collisional properties are here expressed in scattering-relaxation form, making use of the local equilibrium distri¬ 
bution // 9 (x, t). The LB approach to compute the dynamics associated with Eq. |lj) uses sequential computational 
steps. One initially performs a displacement (free-streaming) of each distribution component fi{x) towards the 
nearest-neighbor lattice site pointed at by the discrete velocity u). From there, the equilibrium distribution function 
f^ q (x,t) is computed and the outcome of the collisional process is retrieved. Further iterations of these calculations 
allow the propagation of the lattice dynamics in time. We address the question of whether all these steps can be 
performed in a quantum simulator with practical quantum computing protocols. 

The formal analogy between the Dirac and LB equations was first highlighted in mm » where the velocity dis¬ 
tribution of the particle is treated in a similar fashion as a relativistic spinor. This analogy is further exploited in 
the Majorana representation of the Dirac equation, by using real spinors ]24] . The Dirac (Majorana) equation reads 
(h = 1 here and in the following) 

id t 4b + E,-, (2) 

where we have defined the Dirac (Majorana) streaming matrices qA, mass term /3j j, and the imaginary prefactor i 
proper of quantum mechanical evolution. 
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FIG. 1. (color online) (a) The distribution of the fluid density on a 2-dimensional lattice can be simulated, for example, via 
normal motional modes and internal levels of a set of trapped ions (b). (c) Superposition of two motional modes entangled 
with pseudo spin states can encode velocity distributions in different lattice directions. 


Notice that the streaming matrices of the LB equation are diagonal, while the which generate a Clifford 
algebra, cannot be simultaneously diagonalized. Additionally, the mass matrix (3^ is Hermitian, while standard 
collision matrices come in real symmetric form in the LB equation. Therefore, a complete codification of the LB 
scheme in quantum language requires the implementation of diagonal streaming matrices and of purely imaginary 
symmetric scattering matrices. 

The components of the fluid density distribution function fi(x,t) can be encoded in a set of quantum states 
defined on a proper Fock space. For example, in two dimensions, the distribution of the fluid density over the 
two coordinates can be described by a real quantum wavefunction that encodes the state of two bosonic modes, as 
depicted in Fig. [TJ In the ^-quadrature representation, it reads |\&j) = / dx\dx 2 fi(xi,X 2 ) |xi) \x 2 ), where /*( x±,X 2 ) is 
a real distribution and |a^i( 2 )) the eigenstate of the quadrature of the first (second) bosonic mode. Several quantum 
distributions |\Eh) can be used by entangling the bosonic state to a multi-level system, such as a set of pseudospins, 
therefore the state of the complete system is given by |\I/) = JA 7 ?, |i) <g> (Th), with 77 * being real-valued coefficients. In 
order to keep a real-valued representation of |'F), to be identified with a fluid density distribution function, one has 
to act only with purely imaginary interaction matrices. 

The quantum simulation of the Dirac equation was originally proposed I25| and afterwards realized in a trapped-ion 
experiment [22]. In general, streaming interactions involving matrices in the Dirac or Majorana representation qA V b 
can be implemented by using a pair of pseudospins coupled to one or more bosonic modes. In terms of creation and 
annihilation operators a b (al ) for the bosonic mode in the b direction, one can then consider p b = iS7 b = i{a b — a' b ) 
and write Eq. § on the pseudospin-bosonic Hilbert space of |'F), 

id t (i)> = K b a b i(a b - 4) |*(i)> + 0 l'F(f)}, (3) 

where K b stands for the pseudospin-boson coupling and a b act upon the pseudospin degrees of freedom. 

Thus, the three streaming matrices a b j are written in the Dirac representation as a b = —erf G> a in a pseudospin 
representation and the diagonal mass term as = afI- 2 - These streaming matrices are diagonalized via the operators 
S b = l/y/2(0 + a b ) [53], which have to be physically implemented as quantum gates. Defining S b = exp (—iH b t), 
the associated generators read H b = Aa\ ® I 2 + Baf ® a \, with A = and B = In this way, a purely 

imaginary streaming step i/3V b can be built, which mimics the diagonal streaming of the LB equation. The total 
wavefunction after the streaming steps can be retrieved with a sequential implementation, following the operator 
splitting method [23| . For example, in a 2-dimensional lattice, one has 

l*(tn +1 )> = (S^DySyXS-'DxSJCMtn)) . 


(4) 
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FIG. 2. (Color online) Probabilities of success P for the implementation of the scattering process, (a) Probability of success P a 
per time step of simulating real symmetric random matrices as a function of the number of ancillary measurements N (solid 
lines), together with accumulated probabilities for the whole protocol (dashed lines). Each curve represents a different 
instance of a random matrix, (b) Probability of success of a single step as a function of 7 / 70 , when N = 10. 


The last collision step C , which scrambles particle distributions in different directions, is discussed below. 

Standard collision operators in LB theory are represented by real symmetric matrices associated with non-unitary 
evolution operators. On the other hand, typical controlled quantum mechanics experiments produce unitary dynamics. 
Nevertheless, one can probabilistically encode non-unitary dynamics in a quantum device with a heralded protocol, 
by performing controlled operations conditioned on the state of an ancillary qubit, and then using the state of the 
latter as a flag for the success of the protocol. We consider a purely imaginary symmetric scattering matrix Q, whose 
quantum evolution equation reads idt^i = fproviding a non-unitary evolution operator that describes lattice 
collisions C = exp(— iSlAt). 

The collision operator can be decomposed in a weighted sum of two commuting unitary operators, C = U a + 7 C/ 3 , 
with the constraint 1 1(7(1 < 1 + 7 , assuming without loss of generality that 7 > 0 . 

Given a specific diagonalizable collision operator C and weight 7 , one can then find its decomposition in terms of 
unitaries. In order to find a decomposition in terms of unitaries, C must first be diagonalized as C = VDV'. This 
reduces the problem of finding U a and Up down to an eigenvalue equation, Si = a* + 7 A 7 with Si , a* and /3* being the 
*th eigenvalues of the collision and unitary operators respectively. Notice that, due to the properties of the scattering 
matrix, Si £ R + . Taking into account the normalization conditions, one has the system of equations 


{ Si = at + 7 Pi 

|ai| = l (5) 


The eigenvalues 07 , /?* can now be written as a function of the initial collision operator and weight 7 , 


Re(tti) 


Im(aj) 
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( 6 ) 


The unitary 
the range of 


operators U a ( f 3 ) are reconstructed via = V^ n a n {!3 n )V n j. The 

validity of the method developed here. Simple algebra leads to the set 


real domain of Eqs. ( 6 ) provides 
of inequalities 


| — 1 + <5*| < 7 7 1 + <5*, VI. 


(7) 


By defining 8m and 5 m as the maximal and minimal eigenvalues of the spectrum of (7, the system of inequalities 
in Eq. ([ 7 ]) can be reduced to one of the two inequalities | — 1 + <5 m | < 7 < 1 + <5 m or | — 1 + <5m| < 7 < l + <5 m , 
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respectively when | — 1 + S m \ < \ — 1 + Sm\ or | — 1 + 5 m \ > | — 1 + 5m\- If longer evolution times t are considered, 
the spectral range of C changes accordingly. The weighted 7 -sum derived here can be implemented with quantum 
computing algorithms, using ancillary qubits and controlled U a and Up gates m- By measuring the ancilla state, 
one can determine whether the desired operation has been performed or not. The success of the protocol depends on 
the weighted sum of unitary operators, with a failure probability Pf = y|| U a — U t 3|| 2 /(7 + l) 2 - 

As Pf is an increasing function of 7 , choosing 70 = min{ | — 1 + 6 m |, | — 1 + 5m \} maximizes the probability of success. 
This directly connects the simulation time of the scattering process C with the best choice for 7 . To propagate the 
dynamics of a given collision process C, one can split the step time At into N time intervals A t/N and perform 
the heralded protocol at each step, such that C = exp(—*f2y At/N) N . At each step, one has a collision operator 
exp(—ifiy At/N), with an optimal 70 . In this way, as the step size gets smaller, the success probabilities for each step 
increase, while the total success probability accumulates single success rates from the individual steps. In Fig. [2 |j,, we 
plot the success probability P S (N) = 1 — Pf(N) of the simulation of the single step, as a function of N, for random 
symmetric purely imaginary matrices. As expected, the success probability per step increases as the size for the single 
time step gets smaller. The success of the whole protocol P^ is constant and does not depend on N. In Fig. is 
shown that the optimal protocol is performed at 7 = j 0 . 


DISCUSSION 

The scheme proposed can be adapted to a variety of transport fluid problems. As an example, we consider the 
implementation of an advection-diffusion process in two spatial dimensions. The dynamics of the transported species, 
e.g. pollutants or bacteria, is described by the equation 

dtp + V- ( pU ) = DAp, (8) 

where p = ^/) 4 _ 1 /* is the scalar field transported by a fluid with space-dependent velocity U = ( U x , U y ) and constant 
diffusivity D. 

The problem in Eq. ([ 8 ]) can be recast in LB form, as in Eq. 0 . The corresponding equilibrium distribution function 
is defined as 


n q =w 


pU ■ Ci 


(9) 


with Wi = 1/4, c 2 = 1/2. Note that, by definition, the space-time dependence of the local equilibria is entirely carried 
by the macroscopic fields p and U. 

The scattering matrix reads Ay = Y^t=i A-^WfcA^, where A- 1 ^ = L; = (1,1,1,1), A.- 2 ^ = d x = (1,0,—1,0), 

A- 3 ^ = c iy = ( 0 , 1 , 0 , — 1 ) and A- 4 ^ = c 2 x — cf = ( 1 / 2 , 0 , 1 / 2 , 0 ) are the four eigenvectors. 

The first three corresponding eigenvalues are given by 


Wi =0, UJ 2 = co 3 = 


1 

1/2 + PI cl ’ 


( 10 ) 


which follows from mass conservation and the expression of the diffusion constant D = c 2 s {1/uj — 1/2), respectively. 
By choosing different values for W 2 and W 3 , one can implement anisotropic diffusivities along the x and y directions. 
Finally, w 4 dictates the decay of higher order fields of no direct macroscopic significance, hence it is chosen as W 4 = 1 
so as to annihilate the corresponding held in a single collision step. 

The relative scattering matrix Qy is defined by —Ay [//(of; t) — t)] = Uy/j(i, t). By choosing a Couette how, 

e.g. U = Uo{y, 0), where L is the typical size of the fluid domain, one has = uiip( 1 + Ui), with u\ = — U 3 = U$y/c 2 s 
and U 2 = U 4 = 0. Here, velocities are numbered 1-^4 counterclockwise starting from the +x direction. 

The latter dehnes the quantum scattering matrix as composed of three contributions, namely iffy// = — zAy [fj + 
Wjp 4- Wjpuj\, where iq = — U 3 = Uo/c 2 (a y + a^ y ) is proportional to the position quadrature of the bosonic mode 
associated with the y direction. The three contributions to the scattering matrix represent classical linear wave 
propagation and damping, mass conservation and macroscopic advection, respectively. They can be implemented 
with the quantum simulation protocol previously introduced. The bounds to 7 can be obtained, e.g., for the hrst 
contribution to the scattering matrix —Ay, by computing the spectrum of C = e~ AAt for different time steps At, for 
D = 0.05, in units of I/ 0 C 4 . The result is shown in Fig. [ 3 ] 

Natural quantum platforms for prospective implementation of the proposed scheme could be ions trapped in linear 
Paul traps or superconducting circuit setups, in which the sequential streaming and collision steps in Eq. Q can 
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FIG. 3. (Color online) Spectrum of a collision operator (solid red line) for advection-diffusion process of a four-speed lattice 
as a function of the evolution time step At, in units of 1 / 014 . The allowed region for 7 is bounded by dashed blue lines using 
Eq. Q and shadowed in the picture. 


be realized. The pseudospin-bosonic state can be encoded, in the case of ion traps, in the internal level and motion 
modes of the ions [21], while in a superconducting architectures, one can use the first levels of charge-like qubits, 
e.g. transmon qubits, and microwave resonators [29] . One may consider opening similar avenues in other quantum 
technologies as is the case of quantum photonics PH] and Bose-Einstein condensates [5Tj . 

A practical implementation of the protocol proposed can make use of many-body interactions, involving couplings 
with bosonic modes. These type of gates have been considered in superconducting architectures [321 or in ion-trap 
platforms [33] . For a four-speed lattice, the diagonal streaming processes can be realized with a combination of 
a qubit-boson interaction and two entangling gates among the qubits. For example, the corresponding evolution 
operator for the streaming in the X direction can be written as 


U x = exp 


0a 1 (ai - 4)1 = -Ri t (V 4 )- R 2 t (- 7r / 4 ) c/ c( 7r / 4 ) ex P [K(«i ^ 4)] Uc (tt/ 4 ) R* (tt/ 4 )( —7r / 4 ), ( n ) 


where we have defined an entangling operation between the two qubits Uq = exp [—z(7r/4)(Tj : (j|] and local rotations 
of the i-tli qubit about the j-th axis, 1?/ (0) = expThe U x interaction can then be diagonalized in the 
qubit space via the realization of two S x matrix, S' X U X S X , which can be achieved by a combination of entangling 
two-qubit gates and a phase gate. In the case of more internal degrees of freedom of the lattice, the two-body 


-*(tt / 4 )Ei 


<3 ^ a o 


Similar reasoning 


entangling gate can be substituted by a collective interaction Uc —► exp 
applies to the streaming in the Y direction, considering a different bosonic mode a y . The unitary matrices that 
implement the collision process U a ^ = expcan be implemented in a controlled way [27] by using an 
additional ancillary qubit $.4 and performing the quantum gate exp [ — ^ 1a) ® H a (b)t\- These gates can be 
decomposed in general with a Lie-Trotter-Suzuki decomposition, in terms of many-body interactions of the type 
0 / 4 “of* • • • (a z A ±t A )®H a ( b) = °i“ " -v k N 1 with {ia,ja—k a } £ {x,y, zj. These many-body interactions 

can be obtained by sequential implementation of collective gates and single qubit rotations j32j [33]. In the case of a 
Couette flow, the term with the linear spatial dependence of the scattering matrix can be implemented by considering 
a single qubit rotation entangled with a bosonic displacement, similar to Eq. (11). The quantum resources necessary 
to implement Lie-Trotter-Suzuki decompositions scale polynomially in the internal degrees of freedom of the lattice, 
and sub-polynomially in the digital error [34] . Notice also that the quantum resources needed are invariant with 
respect to the size of the simulated lattice. The latter will depend on the accessibility and readability of distributions 
over Fock spaces in practical implementations, e.g. the ability to characterize distributions over current quadrature 
in superconducting architectures [351 5B . 

Note that the above scheme readily extends to the case of reactive flow, by augmenting the collision operator with 
a local source term proportional to the chemical reaction rate. Such kind of advection-diffusion-reaction phenomena 
in complex geometries, say catalytic reactors, represent a very active area of applications of the LB scheme. Further 
developments may include the implementation of hydrodynamic non-linearities to model the Navier-Stokes fluid 
dynamic equations. This requires the inclusion of quadratic terms in the LB equilibrium distribution. Such nonlinear 
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behavior can be provided in a quantum mechanical experiment by preparing multiple copies of the system mi 
feedback mechanisms [38], or non-unitary operations induced by measurements. 

We have developed a protocol to reproduce the dynamics of fluid transport phenomena in a quantum mechanical 
experiment, by using pseudospins coupled to bosonic modes that can be implemented in different quantum platforms. 
This proposal paves the way to quantum simulation and retrieval of complex classical fluid dynamics in controlled 
quantum systems. 
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